function [Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn,iterVF]=solve_value_function_A3_C(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n)


r=0.036;

% Clenshaw-Curtis quadrature
% n+1 nodes: x in [-1,1] by decreasing order
 [x,wcc] = ccquad(n-1); x=flipud(x);

transition_rates=Par1(1:16);
theta=Par1(21);
a=Par1(22);
b1 = Par1(23); b2 = Par1(24);
gamma=0; 
 

df_1=transition_rates(1);
di_1=transition_rates(2);
lnf_1=transition_rates(3);
lni_1=transition_rates(4);
lfi_1=transition_rates(5);
lif_1=transition_rates(6);
df_2=transition_rates(7);
di_2=transition_rates(8);
lnf_2=transition_rates(9);
lni_2=transition_rates(10);
lfi_2=transition_rates(11);
lif_2=transition_rates(12);
p_1=transition_rates(13);
q_1=transition_rates(14);
p_2=transition_rates(15);
q_2=transition_rates(16);

%preallocation to speed up

derWff = zeros(n); derWfi = zeros(n);derWif = zeros(n);derWii = zeros(n);
derWfn = zeros(1,n);derWin = zeros(1,n);derWnf = zeros(1,n);derWni = zeros(1,n);
dWff_Wfn = zeros(n,1); dWff_Wif = zeros(n,1);dWff_Wnf = zeros(n,1);
dWfi_Wfn = zeros(n,1);dWfi_Wii = zeros(n,1);dWfi_Wnf = zeros(n,1);
dWfi_Wni = zeros(n,1);dWfn_Wif = zeros(n,1);dWfn_Wii = zeros(n,1);
dWfn_Wnf = zeros(n,1);dWfn_Wni = zeros(n,1);dWfn_Wnn = zeros(n,1);
dWif_Wff = zeros(n,1);dWif_Wfn = zeros(n,1);dWif_Wnf = zeros(n,1);
dWii_Wfi = zeros(n,1);dWii_Wfn = zeros(n,1);dWii_Wnf = zeros(n,1);
dWii_Wni = zeros(n,1);dWin_Wfn = zeros(n,1);dWin_Wnf = zeros(n,1);
dWin_Wni = zeros(n,1);dWin_Wnn = zeros(n,1);dWnf_Wfi = zeros(n,1);
dWnf_Wfn = zeros(n,1);dWnf_Wii = zeros(n,1);dWnf_Wni = zeros(n,1);
dWnf_Wnn = zeros(n,1);dWni_Wfn = zeros(n,1);dWni_Wnf = zeros(n,1);
dWni_Wnn = zeros(n,1);
dWff_Wfi = zeros(1,n);dWfi_Wff = zeros(1,n);dWif_Wii = zeros(1,n);
dWif_Win = zeros(1,n);dWii_Wif = zeros(1,n);dWii_Win = zeros(1,n);
eWff = zeros(n);eWfi = zeros(n);eWif = zeros(n);eWii = zeros(n);
eWfn = zeros(1,n);eWin = zeros(1,n);eWnf = zeros(1,n);eWni = zeros(1,n);
 Ff_1_minRw_ff_if_fn_if = zeros(n); Ff_1_minRw_ff_nf_fn_nf = zeros(1,n);
Ff_1_minRw_fi_ii_fn_ii = zeros(n); Ff_1_minRw_fi_ni_fn_ni = zeros(1,n);
Ff_1_Rw_ff_fn = zeros(1,n);Ff_1_Rw_ff_nf = zeros(1,n);Ff_1_Rw_fi_fn = zeros(1,n);
Ff_1_Rw_fi_ni = zeros(1,n);Ff_1_Rw_fn_in = zeros(1,n);Ff_1_Rw_fn_nf = zeros(1,n);
Ff_1_Rw_ff_if = zeros(n); Ff_1_Rw_fi_ii = zeros(n); Ff_1_Rw_fn_nn = zeros(1,n);
Ff_1_Rw_fn_if = zeros(n);Ff_1_Rw_fn_ii = zeros(n);
Ff_2_minRw_ff_fi_nf_fi = zeros(n); Ff_2_minRw_ff_fn_nf_fn = zeros(1,n);
Ff_2_minRw_if_ii_nf_ii = zeros(n); Ff_2_minRw_if_in_nf_in = zeros(1,n);
Ff_2_Rw_ff_nf = zeros(1,n);Ff_2_Rw_ff_fn = zeros(1,n);Ff_2_Rw_if_nf = zeros(1,n);
Ff_2_Rw_if_in = zeros(1,n);Ff_2_Rw_nf_ni = zeros(1,n);Ff_2_Rw_nf_fn = zeros(1,n);
Ff_2_Rw_ff_fi = zeros(n); Ff_2_Rw_if_ii = zeros(n); Ff_2_Rw_nf_in = zeros(1,n);
Ff_2_Rw_nf_fi = zeros(n);Ff_2_Rw_nf_ii = zeros(n);
Fi_1_Rw_if_ff = zeros(n);Fi_1_Rw_ii_fi = zeros(n);
Fi_2_Rw_fi_ff = zeros(n);Fi_1_Rw_ii_if = zeros(n);
Fi_1_Rw_if_nf = zeros(1,n); Fi_1_Rw_ii_ni = zeros(1,n);
Fi_1_Rw_in_fn = zeros(1,n); Fi_1_Rw_in_nf = zeros(1,n);
Fi_1_Rw_in_ni = zeros(1,n); Fi_2_Rw_fi_fn = zeros(1,n);
Fi_1_Rw_ii_in = zeros(1,n); Fi_1_Rw_ni_fn = zeros(1,n);
Fi_1_Rw_ni_in = zeros(1,n); Fi_1_Rw_ni_nf = zeros(1,n);
intmaxWff_Wfn_Wif = zeros(n); intmaxWff_Wnf_Wfi = zeros(n);
intmaxWfi_Wff = zeros(n);intmaxWfi_Wfn_Wii = zeros(n);intmaxWif_Wff = zeros(n);
intmaxWif_Wnf_Wii = zeros(n);intmaxWii_Wfi = zeros(n);intmaxWii_Wif = zeros(n);
intmaxWff_Wfn_Wnf = zeros(1,n);intmaxWff_Wnf_Wfn = zeros(1,n);
intmaxWfi_Wfn = zeros(1,n);intmaxWfi_Wfn_Wni = zeros(1,n);
intmaxWif_Wnf_Win = zeros(1,n);intmaxWii_Win = zeros(1,n);
intmaxWfn_Win = zeros(1,n);intmaxWif_Wnf = zeros(1,n);
intmaxWii_Wni = zeros(1,n);intmaxWin_Wfn = zeros(1,n);
intmaxWin_Wnf = zeros(1,n);intmaxWin_Wni = zeros(1,n);
intmaxWnf_Wni = zeros(1,n);intmaxWni_Wfn = zeros(1,n);
intmaxWni_Win = zeros(1,n);intmaxWni_Wnf = zeros(1,n);
minRw_fi_ff_fi_nf = zeros(n); minRw_if_ff_if_fn = zeros(n);
minRw_ii_fi_ii_fn = zeros(n);minRw_ii_if_ii_nf = zeros(n);
Rw_ff_fi = zeros(n);Rw_ff_if = zeros(n);Rw_fi_ff = zeros(n);Rw_fi_ii = zeros(n);
Rw_fn_if = zeros(n);Rw_fn_ii = zeros(n);Rw_if_ff = zeros(n);Rw_if_ii = zeros(n);
Rw_ii_fi = zeros(n);Rw_ii_if = zeros(n);Rw_nf_fi = zeros(n);Rw_nf_ii = zeros(n);
Rw_ff_fn = zeros(1,n);Rw_ff_nf = zeros(1,n);Rw_fi_fn = zeros(1,n);Rw_fi_ni = zeros(1,n);
Rw_fn_in = zeros(1,n);Rw_fn_nf = zeros(1,n);Rw_fn_ni = zeros(1,n);Rw_if_in = zeros(1,n);
Rw_if_nf = zeros(1,n);Rw_ii_in = zeros(1,n);Rw_ii_ni = zeros(1,n);Rw_in_fn = zeros(1,n);
Rw_in_nf = zeros(1,n);Rw_in_ni = zeros(1,n);Rw_nf_fn = zeros(1,n);Rw_nf_in = zeros(1,n);
Rw_nf_ni = zeros(1,n);Rw_ni_fn = zeros(1,n);Rw_ni_in = zeros(1,n);Rw_ni_nf = zeros(1,n);
 Wff = zeros(n);
Wff_max = zeros(n,1); Wff_min = zeros(n,1); 
Wff_Rw_ff_fi = zeros(n);Wff_Rw_ff_fn = zeros(1,n);
Wff_Rw_ff_if = zeros(n);Wff_Rw_ff_nf = zeros(1,n);
Wfi = zeros(n);
Wfi_max = zeros(n,1); Wfi_min = zeros(n,1); 
Wfi_Rw_fi_ff = zeros(n);Wfi_Rw_fi_fn = zeros(1,n);
Wfi_Rw_fi_ii = zeros(n);Wfi_Rw_fi_ni = zeros(1,n);
Wfn = zeros(n,1); Wfn_Rw_fn_if = zeros(n); Wfn_Rw_fn_ii = zeros(n);
Wfn_Rw_fn_in = zeros(1,n);Wfn_Rw_fn_nf = zeros(1,n);Wfn_Rw_fn_ni = zeros(1,n);
Wif = zeros(n);
Wif_max = zeros(n,1); Wif_min = zeros(n,1); Wif_Rw_if_ff = zeros(n);
Wif_Rw_if_ii = zeros(n);Wif_Rw_if_in = zeros(1,n);Wif_Rw_if_nf = zeros(1,n);
Wii = zeros(n);
Wii_max = zeros(n,1); Wii_min = zeros(n,1); Wii_Rw_ii_fi = zeros(n);
Wii_Rw_ii_if = zeros(n);Wii_Rw_ii_in = zeros(1,n);Wif_Rw_ii_ni = zeros(1,n);
Win = zeros(n,1); Win_Rw_in_fn=zeros(1,n); Win_Rw_in_nf=zeros(1,n); Win_Rw_in_ni=zeros(1,n);
Wnf = zeros(n,1); Wnf_Rw_nf_fn=zeros(1,n); Wnf_Rw_nf_in=zeros(1,n); Wnf_Rw_nf_ni=zeros(1,n);
Wnf_Rw_nf_fi=zeros(n); Wnf_Rw_nf_ii=zeros(n);
Wni = zeros(n,1); Wni_Rw_ni_fn=zeros(1,n); Win_Rw_ni_in=zeros(1,n); Win_Rw_ni_nf=zeros(1,n);
%% Transforming the probabilities g and f in densities
% If wages are linearly spaced
wht_f_1=(wf_1(n)-wf_1(1))/(n-1); 
wht_i_1=(wi_1(n)-wi_1(1))/(n-1);
wht_f_2=(wf_2(n)-wf_2(1))/(n-1);
wht_i_2=(wi_2(n)-wi_2(1))/(n-1);

    % Initializing
    
    % value function (min and max)
    Wnn=(1-exp(-theta*(min(wi_1)+min(wi_2))))/theta/r;
    Wfn_min= (1-exp(-theta*(min(wf_1)+min(wi_2))))/theta/r; Wfn_max= (1-exp(-theta*(max(wf_1)+max(wi_2))))/theta/r;
    Win_min=(1-exp(-theta*(min(wi_1)+min(wi_2))))/theta/r; Win_max=(1-exp(-theta*(max(wi_1)+max(wi_2))))/theta/r ;
    Wnf_min=(1-exp(-theta*(min(wi_1)+min(wf_2))))/theta/r ; Wnf_max=(1-exp(-theta*(max(wi_1)+max(wf_2))))/theta/r ;
    Wni_min=(1-exp(-theta*(min(wi_1)+min(wi_2))))/theta/r ; Wni_max=(1-exp(-theta*(max(wi_1)+max(wi_2))))/theta/r ;
    Wff_min=(1-exp(-theta*(min(wf_1)+wf_2)))/theta/r; Wff_max= (1-exp(-theta*(max(wf_1)+wf_2)))/theta/r;
    Wii_min=(1-exp(-theta*(min(wi_1)+wi_2)))/theta/r ; Wii_max= (1-exp(-theta*(max(wi_1)+wi_2)))/theta/r ;
    Wfi_min=(1-exp(-theta*(min(wf_1)+wi_2)))/theta/r ; Wfi_max=(1-exp(-theta*(max(wf_1)+wi_2)))/theta/r;
    Wif_min=(1-exp(-theta*(min(wi_1)+wf_2)))/theta/r ; Wif_max=(1-exp(-theta*(max(wi_1)+wf_2)))/theta/r;
    
    Wfn=(Wfn_min+Wfn_max)/2+x*(Wfn_max-Wfn_min)/2; % Interpolation with CC quadrature (defined above)
    Win=(Win_min+Win_max)/2+x*(Win_max-Win_min)/2;
    Wnf=(Wnf_min+Wnf_max)/2+x*(Wnf_max-Wnf_min)/2;
    Wni=(Wni_min+Wni_max)/2+x*(Wni_max-Wni_min)/2;
    Wff=(Wff_min+Wff_max)*ones(1,n)/2+(Wff_max-Wff_min)*x'/2;
    Wii=(Wii_min+Wii_max)*ones(1,n)/2+(Wii_max-Wii_min)*x'/2;
    Wfi=(Wfi_min+Wfi_max)*ones(1,n)/2+(Wfi_max-Wfi_min)*x'/2;
    Wif=(Wif_min+Wif_max)*ones(1,n)/2+(Wif_max-Wif_min)*x'/2;
    
    critVF=1;
    iterVF=1;
    omegaW=0.1;
    while mean(critVF)>0.005 && iterVF<=400
        
        % Reservation wages
        for i=1:n
            % to construct Wfn
            dWni_Wfn=abs(Wni-Wfn(i))./Wfn(i);
            Rw_ni_fn(i)=wi_2(min(find(dWni_Wfn<=min(dWni_Wfn))));
            Wni_Rw_ni_fn(i)=Wni(min(find(dWni_Wfn<=min(dWni_Wfn))));
            Fi_2_Rw_ni_fn(i)=Fi_2(min(find(dWni_Wfn<=min(dWni_Wfn))));
            dWin_Wfn=abs(Win-Wfn(i))./Wfn(i);
            Rw_in_fn(i)=wi_1(min(find(dWin_Wfn<=min(dWin_Wfn))));
            Win_Rw_in_fn(i)=Win(min(find(dWin_Wfn<=min(dWin_Wfn))));
            Fi_1_Rw_in_fn(i)=Fi_1(min(find(dWin_Wfn<=min(dWin_Wfn))));
            dWff_Wfn=abs(Wff(i,:)-Wfn(i))./Wfn(i);
            Rw_ff_fn(i)=wf_2(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
            Wff_Rw_ff_fn(i)=Wff(i,find(dWff_Wfn<=min(dWff_Wfn), 1 ));
            Ff_2_Rw_ff_fn(i)=Ff_2(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
            dWnf_Wfn=abs(Wnf-Wfn(i))./Wfn(i);
            Rw_nf_fn(i)=wf_2(find(dWnf_Wfn<=min(dWnf_Wfn), 1 ));
            Wnf_Rw_nf_fn(i)=Wnf(min(find(dWnf_Wfn<=min(dWnf_Wfn))));
            Ff_2_Rw_nf_fn(i)=Ff_2(find(dWnf_Wfn<=min(dWnf_Wfn), 1 ));
            Ff_2_minRw_ff_fn_nf_fn(i)=Ff_2(min(find(dWnf_Wfn<=min(dWnf_Wfn), 1 ),min(find(dWff_Wfn<=min(dWff_Wfn)))));
            dWff_Wnf=abs(Wff(i,:)'-Wnf)./Wnf;
            Rw_ff_nf(i)=wf_2(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
            Wff_Rw_ff_nf(i)=Wff(i,find(dWff_Wnf<=min(dWff_Wnf), 1 ));
            Ff_2_Rw_ff_nf(i)=Ff_2(min(find(dWff_Wnf<=min(dWff_Wnf))));
            dWfi_Wfn=abs(Wfi(i,:)-Wfn(i))./Wfn(i);
            Rw_fi_fn(i)=wi_2(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
            Wfi_Rw_fi_fn(i)=Wfi(i,find(dWfi_Wfn<=min(dWfi_Wfn), 1 ));
            Fi_2_Rw_fi_fn(i)=Fi_2(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
            
            derWfn(i)=1./(r+df_1*(1-p_2)+df_1*p_2*Fi_2_Rw_ni_fn(i)+lfi_1*Fi_1_Rw_in_fn(i) ...
                +lnf_2*Ff_2_minRw_ff_fn_nf_fn(i)+lni_2*Fi_2_Rw_fi_fn(i));
            
            % to construct Win
            dWni_Win=abs(Wni(i,:)-Win(i))./Win(i);
            Rw_ni_in(i)=wi_2(min(find(dWni_Win<=min(dWni_Win))));
            Wni_Rw_ni_in(i)=Wni(find(dWni_Win<=min(dWni_Win), 1 ));
            Fi_2_Rw_ni_in(i)=Fi_2(find(dWni_Win<=min(dWni_Win), 1 ));
            dWfn_Win=abs(Wfn(i,:)-Win(i))./Win(i);
            Rw_fn_in(i)=wf_1(find(dWfn_Win<=min(dWfn_Win), 1 ));
            Wfn_Rw_fn_in(i)=Wfn(find(dWfn_Win<=min(dWfn_Win), 1 ));
            Ff_1_Rw_fn_in(i)=Ff_1(find(dWfn_Win<=min(dWfn_Win), 1 ));
            dWif_Win=abs(Wif(i,:)-Win(i))./Win(i);
            Rw_if_in(i)=wf_2(find(dWif_Win<=min(dWif_Win), 1 ));
            Wif_Rw_if_in(i)=Wif(i,find(dWif_Win<=min(dWif_Win), 1 ));
            Ff_2_Rw_if_in(i)=Ff_2(find(dWif_Win<=min(dWif_Win), 1 ));
            dWnf_Win=abs(Wnf(i,:)-Win(i))./Win(i);
            Rw_nf_in(i)=wf_2(find(dWnf_Win<=min(dWnf_Win), 1 ));
            Wnf_Rw_nf_in(i)=Wnf(find(dWnf_Win<=min(dWnf_Win), 1 ));
            Ff_2_Rw_nf_in(i)=Ff_2(find(dWnf_Win<=min(dWnf_Win), 1 ));
            Ff_2_minRw_if_in_nf_in(i)=Ff_2(min(find(dWif_Win<=min(dWif_Win), 1 ),min(find(dWnf_Win<=min(dWnf_Win)))));
            dWif_Wnf=abs(Wif(i,:)'-Wnf)./Wnf;
            Rw_if_nf(i)=wf_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
            Wif_Rw_if_nf(i)=Wif(i,min(find(dWif_Wnf<=min(dWif_Wnf))));
            Ff_2_Rw_if_nf(i)=Ff_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
            dWii_Win=abs(Wii(i,:)-Win(i))./Win(i);
            Rw_ii_in(i)=wi_2(min(find(dWii_Win<=min(dWii_Win))));
            Wii_Rw_ii_in(i)=Wii(i,find(dWii_Win<=min(dWii_Win), 1 ));
            Fi_2_Rw_ii_in(i)=Fi_2(min(find(dWii_Win<=min(dWii_Win))));
            
            
            derWin(i)=1./(r+di_1*(1-q_2)+di_1*q_2*Fi_2_Rw_ni_in(i)+lif_1*Ff_1_Rw_fn_in(i) ...
                +lnf_2*Ff_2_minRw_if_in_nf_in(i)+lni_2*Fi_2_Rw_ii_in(i));
        end
        
        for j=1:n
            % to construct Wnf
            dWin_Wnf=abs(Win-Wnf(j))./Wnf(j);
            Rw_in_nf(j)=wi_1(min(find(dWin_Wnf<=min(dWin_Wnf))));
            Win_Rw_in_nf(j)=Win(min(find(dWin_Wnf<=min(dWin_Wnf))));
            Fi_1_Rw_in_nf(j)=Fi_1(min(find(dWin_Wnf<=min(dWin_Wnf))));
            dWni_Wnf=abs(Wni-Wnf(j))./Wnf(j);
            Rw_ni_nf(j)=wi_2(min(find(dWni_Wnf<=min(dWni_Wnf))));
            Wni_Rw_ni_nf(j)=Wni(min(find(dWni_Wnf<=min(dWni_Wnf))));
            Fi_2_Rw_ni_nf(j)=Fi_2(min(find(dWni_Wnf<=min(dWni_Wnf))));
            dWff_Wnf=abs(Wff(:,j)-Wnf(j))./Wnf(j);
            Rw_ff_nf(j)=wf_1(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
            Wff_Rw_ff_nf(j)=Wff(min(find(dWff_Wnf<=min(dWff_Wnf))),j);
            Ff_1_Rw_ff_nf(j)=Ff_1(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
            dWfn_Wnf=abs(Wfn-Wnf(j))./Wnf(j);
            Rw_fn_nf(j)=wf_1(find(dWfn_Wnf<=min(dWfn_Wnf), 1 ));
            Wfn_Rw_fn_nf(j)=Wfn(find(dWfn_Wnf<=min(dWfn_Wnf), 1 ));
            Ff_1_Rw_fn_nf(j)=Ff_1(min(find(dWfn_Wnf<=min(dWfn_Wnf))));
            Ff_1_minRw_ff_nf_fn_nf(j)=Ff_1(min(min(find(dWff_Wnf<=min(dWff_Wnf))),find(dWfn_Wnf<=min(dWfn_Wnf), 1 )));
            dWff_Wfn=abs(Wff(:,j)-Wfn)./Wfn;
            Rw_ff_fn(j)=wf_1(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
            Wff_Rw_ff_fn(j)=Wff(find(dWff_Wfn<=min(dWff_Wfn), 1 ),j);
            Ff_1_Rw_ff_fn(j)=Ff_1(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
            dWif_Wnf=abs(Wif(:,j)-Wnf(j))./Wnf(j);
            Rw_if_nf(j)=wi_1(min(find(dWif_Wnf<=min(dWif_Wnf))));
            Wif_Rw_if_nf(j)=Wif(min(find(dWif_Wnf<=min(dWif_Wnf))),j);
            Fi_1_Rw_if_nf(j)=Fi_1(find(dWif_Wnf<=min(dWif_Wnf), 1 ));
            % to construct Wni
            dWin_Wni=abs(Win-Wni(j))./Wni(j);
            Rw_in_ni(j)=wi_1(find(dWin_Wni<=min(dWin_Wni), 1 ));
            Win_Rw_in_ni(j)=Win(min(find(dWin_Wni<=min(dWin_Wni))));
            Fi_1_Rw_in_ni(j)=Fi_1(find(dWin_Wni<=min(dWin_Wni), 1 ));
            dWnf_Wni=abs(Wnf-Wni(j))./Wni(j);
            Rw_nf_ni(j)=wf_2(min(find(dWnf_Wni<=min(dWnf_Wni))));
            Wnf_Rw_nf_ni(j)=Wnf(min(find(dWnf_Wni<=min(dWnf_Wni))));
            Ff_2_Rw_nf_ni(j)=Ff_2(find(dWnf_Wni<=min(dWnf_Wni), 1 ));
            dWfi_Wni=abs(Wfi(:,j)-Wni(j))./Wni(j);
            Rw_fi_ni(j)=wf_1(min(find(dWfi_Wni<=min(dWfi_Wni))));
            Wfi_Rw_fi_ni(j)=Wfi(min(find(dWfi_Wni<=min(dWfi_Wni))),j);
            Ff_1_Rw_fi_ni(j)=Ff_1(find(dWfi_Wni<=min(dWfi_Wni), 1 ));
            dWfn_Wni=abs(Wfn-Wni(j))./Wni(j);
            Rw_fn_ni(j)=wf_1(find(dWfn_Wni<=min(dWfn_Wni), 1 ));
            Wfn_Rw_fn_ni(j)=Wfn(find(dWfn_Wni<=min(dWfn_Wni), 1 ));
            Ff_1_Rw_fn_ni(j)=Ff_1(find(dWfn_Wni<=min(dWfn_Wni), 1 ));
            Ff_1_minRw_fi_ni_fn_ni(j)=Ff_1(min(min(find(dWfi_Wni<=min(dWfi_Wni))),min(find(dWfn_Wni<=min(dWfn_Wni)))));
            dWfi_Wfn=abs(Wfi(:,j)-Wfn)./Wfn;
            Rw_fi_fn(j)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
            Wfi_Rw_fi_fn(j)=Wfi(min(find(dWfi_Wfn<=min(dWfi_Wfn))),j);
            Ff_1_Rw_fi_fn(j)=Ff_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
            dWii_Wni=abs(Wii(:,j)-Wni(j))./Wni(j);
             Rw_ii_ni(j)=wi_1(min(find(dWii_Wni<=min(dWii_Wni))));
            Wii_Rw_ii_ni(j)=Wii(min(find(dWii_Wni<=min(dWii_Wni))),j);
            Fi_1_Rw_ii_ni(j)=Fi_1(min(find(dWii_Wni<=min(dWii_Wni))));
            
            derWnf(j)=1./(r+df_2*(1-p_1)+df_2*p_1*Fi_1_Rw_in_nf(j)+lfi_2*Fi_2_Rw_ni_nf(j) ...
                +lnf_1*Ff_1_minRw_ff_nf_fn_nf(j)+lni_1*Fi_1_Rw_if_nf(j));
            
            derWni(j)=1./(r+di_2*(1-q_1)+di_2*q_1*Fi_1_Rw_in_ni(j)+lif_2*Ff_2_Rw_nf_ni(j) ...
                +lnf_1*Ff_1_minRw_fi_ni_fn_ni(j)+lni_1*Fi_1_Rw_ii_ni(j));
            
        end
        
        for i=1:n
            for j=1:n
                
                % to construct Wff
                dWif_Wff=abs(Wif(:,j)-Wff(i,j))./Wff(i,j);dWif_Wff=abs(Wif(:,n)-Wff(i,n))./Wff(i,n);
                Rw_if_ff(i,j)=wi_1(min(find(dWif_Wff<=min(dWif_Wff))));
                Wif_Rw_if_ff(i,j)=Wif(min(find(dWif_Wff<=min(dWif_Wff))),j);
                Fi_1_Rw_if_ff(i,j)=Fi_1(min(find(dWif_Wff<=min(dWif_Wff))));
                dWfi_Wff=abs(Wfi(i,:)-Wff(i,j))./Wff(i,j);dWfi_Wff=abs(Wfi(n,:)-Wff(n,j))./Wff(n,j);
                Rw_fi_ff(i,j)=wi_2(min(find(dWfi_Wff<=min(dWfi_Wff))));
                Wfi_Rw_fi_ff(i,j)=Wfi(i,min(find(dWfi_Wff<=min(dWfi_Wff))));
                Fi_2_Rw_fi_ff(i,j)=Fi_2(min(find(dWfi_Wff<=min(dWfi_Wff))));
                
                derWff(i,j)=1./(r+df_1+df_2+lfi_1*Fi_1_Rw_if_ff(i,j) ...
                    +lfi_2*Fi_2_Rw_fi_ff(i,j));
                
                % to construct Wii
                dWfi_Wii=abs(Wfi(:,j)-Wii(i,j))./Wii(i,j);dWfi_Wii=abs(Wfi(:,n)-Wii(i,n))./Wii(i,n);
                Rw_fi_ii(i,j)=wf_1(min(find(dWfi_Wii<=min(dWfi_Wii))));
                Wfi_Rw_fi_ii(i,j)=Wfi(min(find(dWfi_Wii<=min(dWfi_Wii))),j);
                Ff_1_Rw_fi_ii(i,j)=Ff_1(find(dWfi_Wii<=min(dWfi_Wii), 1 ));
                dWfn_Wii=abs(Wfn-Wii(i,j))./Wii(i,j);
                Rw_fn_ii(i,j)=wf_1(min(find(dWfn_Wii<=min(dWfn_Wii))));
                Wfn_Rw_fn_ii(i,j)=Wfn(find(dWfn_Wii<=min(dWfn_Wii), 1 ));
                Ff_1_Rw_fn_ii(i,j)=Ff_1(min(find(dWfn_Wii<=min(dWfn_Wii))));
                Ff_1_minRw_fi_ii_fn_ii(i,j)=Ff_1(min(min(find(dWfi_Wii<=min(dWfi_Wii))),find(dWfn_Wii<=min(dWfn_Wii), 1 )));
                dWfi_Wfn=abs(Wfi(:,j)-Wfn)./Wfn;
                Rw_fi_fn(j)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
                Wfi_Rw_fi_fn(j)=Wfi(min(find(dWfi_Wfn<=min(dWfi_Wfn))),j);
                Ff_1_Rw_fi_fn(j)=Ff_1(find(dWfi_Wfn<=min(dWfi_Wfn), 1 ));
                dWif_Wii=abs(Wif(i,:)-Wii(i,j))./Wii(i,j);dWif_Wii=abs(Wif(n,:)-Wii(n,j))./Wii(n,j);
                Rw_if_ii(i,j)=wf_2(find(dWif_Wii<=min(dWif_Wii), 1 ));
                Wif_Rw_if_ii(i,j)=Wif(i,find(dWif_Wii<=min(dWif_Wii), 1 ));
                Ff_2_Rw_if_ii(i,j)=Ff_2(find(dWif_Wii<=min(dWif_Wii), 1 ));
                dWnf_Wii=abs(Wnf-Wii(i,j))./Wii(i,j);
                Rw_nf_ii(i,j)=wf_2(min(find(dWnf_Wii<=min(dWnf_Wii))));
                Wnf_Rw_nf_ii(i,j)=Wnf(find(dWnf_Wii<=min(dWnf_Wii), 1 ));
                Ff_2_Rw_nf_ii(i,j)=Ff_2(find(dWnf_Wii<=min(dWnf_Wii), 1 ));
                Ff_2_minRw_if_ii_nf_ii(i,j)=Ff_2(min(find(dWif_Wii<=min(dWif_Wii), 1 ),min(find(dWnf_Wii<=min(dWnf_Wii)))));
                dWif_Wnf=abs(Wif(i,:)'-Wnf)./Wnf;
                Rw_if_nf(j)=wf_2(find(dWif_Wnf<=min(dWif_Wnf), 1 ));
                 Wif_Rw_if_nf(j)=Wif(i,find(dWif_Wnf<=min(dWif_Wnf), 1 ));
                Ff_2_Rw_if_nf(j)=Ff_2(find(dWif_Wnf<=min(dWif_Wnf), 1 ));
                
                derWii(i,j)=1./(r+di_1+di_2+lif_1*Ff_1_minRw_fi_ii_fn_ii(i,j) ...
                    +lif_2*Ff_2_minRw_if_ii_nf_ii(i,j));
                
                % to construct Wfi
                dWii_Wfi=abs(Wii(:,j)-Wfi(i,j))./Wfi(i,j);dWii_Wfi=abs(Wii(:,n)-Wfi(i,n))./Wfi(i,n);
                Rw_ii_fi(i,j)=wi_1(min(find(dWii_Wfi<=min(dWii_Wfi))));
                Wii_Rw_ii_fi(i,j)=Wii(min(find(dWii_Wfi<=min(dWii_Wfi))),j);
                Fi_1_Rw_ii_fi(i,j)=Fi_1(find(dWii_Wfi<=min(dWii_Wfi), 1 ));
                dWff_Wfi=abs(Wff(i,:)-Wfi(i,j))./Wfi(i,j);dWff_Wfi=abs(Wff(n,:)-Wfi(n,j))./Wfi(n,j);
                Rw_ff_fi(i,j)=wf_2(min(find(dWff_Wfi<=min(dWff_Wfi))));
                Wff_Rw_ff_fi(i,j)=Wff(i,min(find(dWff_Wfi<=min(dWff_Wfi))));
                Ff_2_Rw_ff_fi(i,j)=Ff_2(min(find(dWff_Wfi<=min(dWff_Wfi))));
                dWnf_Wfi=abs(Wnf-Wfi(i,j))./Wfi(i,j);
                Rw_nf_fi(i,j)=wf_2(find(dWnf_Wfi<=min(dWnf_Wfi), 1 ));
                Wnf_Rw_nf_fi(i,j)=Wnf(find(dWnf_Wfi<=min(dWnf_Wfi), 1 ));
                Ff_2_Rw_nf_fi(i,j)=Ff_2(min(find(dWnf_Wfi<=min(dWnf_Wfi))));
                Ff_2_minRw_ff_fi_nf_fi(i,j)=Ff_2(min(min(find(dWff_Wfi<=min(dWff_Wfi))),min(find(dWnf_Wfi<=min(dWnf_Wfi)))));
                dWff_Wnf=abs(Wff(i,:)'-Wnf)./Wnf;
                Rw_ff_nf(i)=wf_2(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
                Wff_Rw_ff_nf(i)=Wff(i,min(find(dWff_Wnf<=min(dWff_Wnf))));
                Ff_2_Rw_ff_nf(i)=Ff_2(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
                
                derWfi(i,j)=1./(r+df_1+di_2+lfi_1*Fi_1_Rw_ii_fi(i,j) ...
                    +lif_2*Ff_2_minRw_ff_fi_nf_fi(i,j));
                
                % to construct Wif
                dWff_Wif=abs(Wff(:,j)-Wif(i,j))./Wif(i,j);dWff_Wif=abs(Wff(:,n)-Wif(i,n))./Wif(i,n);
                Rw_ff_if(i,j)=wf_1(min(find(dWff_Wif<=min(dWff_Wif))));
                Wff_Rw_ff_if(i,j)=Wff(min(find(dWff_Wif<=min(dWff_Wif))),j);
                Ff_1_Rw_ff_if(i,j)=Ff_1(min(find(dWff_Wif<=min(dWff_Wif))));
                dWfn_Wif=abs(Wfn-Wif(i,j))./Wif(i,j);
                Rw_fn_if(i,j)=wf_1(min(find(dWfn_Wif<=min(dWfn_Wif))));
                Wfn_Rw_fn_if(i,j)=Wfn(min(find(dWfn_Wif<=min(dWfn_Wif))));
                Ff_1_Rw_fn_if(i,j)=Ff_1(find(dWfn_Wif<=min(dWfn_Wif), 1 ));
                Ff_1_minRw_ff_if_fn_if(i,j)=Ff_1(min(find(dWff_Wif<=min(dWff_Wif), 1 ),find(dWfn_Wif<=min(dWfn_Wif), 1 )));
                dWff_Wfn=abs(Wff(:,j)-Wfn)./Wfn;
                Rw_ff_fn(j)=wf_1(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
                Wff_Rw_ff_fn(j)=Wff(find(dWff_Wfn<=min(dWff_Wfn), 1 ),j);
                Ff_1_Rw_ff_fn(j)=Ff_1(min(find(dWff_Wfn<=min(dWff_Wfn))));
                dWii_Wif=abs(Wii(i,:)-Wif(i,j))./Wif(i,j);dWii_Wif=abs(Wii(n,:)-Wif(n,j))./Wif(n,j);
                Rw_ii_if(i,j)=wi_2(min(find(dWii_Wif<=min(dWii_Wif))));
                Wii_Rw_ii_if(i,j)=Wii(min(find(dWii_Wif<=min(dWii_Wif))),j);
                Fi_2_Rw_ii_if(i,j)=Fi_2(find(dWii_Wif<=min(dWii_Wif), 1 ));
                
                derWif(i,j)=1./(r+di_1+df_2+lif_1*Ff_1_minRw_ff_if_fn_if(i,j) ...
                    +lfi_2*Fi_2_Rw_ii_if(i,j));
                
                % to construct Wnn
                dWfn_Wnn=abs(Wfn-Wnn)./Wnn;
                Rw_fn_nn=wf_1(find(dWfn_Wnn<=min(dWfn_Wnn), 1 ));
                Ff_1_Rw_fn_nn=Ff_1(min(find(dWfn_Wnn<=min(dWfn_Wnn))));
                
                dWnf_Wnn=abs(Wnf-Wnn)./Wnn;
                Rw_nf_nn=wf_2(min(find(dWnf_Wnn<=min(dWnf_Wnn))));
                Ff_2_Rw_nf_nn=Ff_2(min(find(dWnf_Wnn<=min(dWnf_Wnn))));
               
                dWin_Wnn=abs(Win-Wnn)./Wnn;
                Rw_in_nn=wi_1(find(dWin_Wnn<=min(dWin_Wnn), 1 ));
                Fi_1_Rw_in_nn=Fi_1(min(find(dWin_Wnn<=min(dWin_Wnn))));
                
                dWni_Wnn=abs(Wni-Wnn)./Wnn;
                Rw_ni_nn=wi_2(find(dWni_Wnn<=min(dWni_Wnn), 1 ));
                xFi_2_Rw_ni_nn=Fi_2(find(dWni_Wnn<=min(dWni_Wnn), 1 ));
                
                % to construct min{wij-fj,wij-fn} - for spouse 1
                 dWii_Wfn=abs(Wii(:,j)-Wfn)./Wfn;
                minRw_ii_fi_ii_fn(i,j)=wi_1(min(min(find(dWii_Wfi<=min(dWii_Wfi))),min(find(dWii_Wfn<=min(dWii_Wfn)))));
                dWif_Wfn=abs(Wif(:,j)-Wfn)./Wfn;
                minRw_if_ff_if_fn(i,j)=wi_1(min(find(dWif_Wff<=min(dWif_Wff), 1 ),find(dWif_Wfn<=min(dWif_Wfn), 1 )));
                
                % to construct min{wji-jf,wji-nf} - for spouse 2
                dWii_Wnf=abs(Wii(i,:)'-Wnf)./Wnf;
                minRw_ii_if_ii_nf(i,j)=wi_2(min(min(find(dWii_Wif<=min(dWii_Wif))),min(find(dWii_Wnf<=min(dWii_Wnf)))));
                dWfi_Wnf=abs(Wfi(i,:)'-Wnf)./Wnf;
                minRw_fi_ff_fi_nf(i,j)=wi_2(min(min(find(dWfi_Wff<=min(dWfi_Wff))),min(find(dWfi_Wnf<=min(dWfi_Wnf)))));
                
            end
        end
        
        % solve integral by parts
        
        % Wfn
        for o=1:n
            dWff_Wfn=abs(Wff(o,:)-Wfn(o))./Wfn(o);
            Rw_ff_fn(o)=wf_2(min(find(dWff_Wfn<=min(dWff_Wfn))));
            dWff_Wnf=abs(Wff(o,:)'-Wnf)./Wnf;
            Rw_ff_nf(o)=wf_2(min(find(dWff_Wnf<=min(dWff_Wnf))));
            dWfi_Wfn=abs(Wfi(o,:)-Wfn(o))./Wfn(o);
            Rw_fi_fn(o)=wi_2(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        end
        for i=1:n
            
            intmaxWni_Wfn(i)=sum((wi_2>Rw_ni_fn(i)).*Fi_2.*derWni'.*wht_i_2);
            intmaxWin_Wfn(i)=sum((wi_1>Rw_in_fn(i)).*Fi_1.*derWin'.*wht_i_1);
            intmaxWfi_Wfn(i)=sum((wi_2>Rw_fi_fn(i)).*Fi_2.*derWfi(i,:)'.*wht_i_2);
            %intmaxWff_Wnf_Wfn(i)=sum((wf_2>Rw_ff_fn(i)).*Ff_2.*derWff(i,:)'.*wht_f_2);
            intmaxWff_Wnf_Wfn(i)=max(0,sum((wf_2<Rw_ff_nf(i)).*Ff_2.*derWff(i,:)'.*wht_f_2+(wf_2>=Rw_ff_nf(i)).*Ff_2.*derWnf'.*wht_f_2));
            
        end
        % Win
        
        for o=1:n
            dWif_Wnf=abs(Wif(o,:)'-Wnf)./Wnf;
            Rw_if_nf(o)=wf_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
        end
        
        for i=1:n
            
            intmaxWni_Win(i)=sum((wi_2>Rw_ni_in(i)).*Fi_2.*derWni'.*wht_i_2);
            intmaxWfn_Win(i)=sum((wf_1>Rw_fn_in(i)).*Ff_1.*derWfn'.*wht_f_1);
            intmaxWii_Win(i)=sum((wi_2>Rw_ii_in(i)).*Fi_2.*derWii(i,:)'.*wht_i_2);
            %intmaxWif_Wnf_Win(i)=sum((wf_2>Rw_if_in(i)).*Ff_2.*derWif(i,:)'.*wht_f_2);
            intmaxWif_Wnf_Win(i)=max(0,sum((wf_2<Rw_if_nf(i)).*Ff_2.*derWif(i,:)'.*wht_f_2+(wf_2>=Rw_if_nf(i)).*Ff_2.*derWnf'.*wht_f_2));
            
        end

        
        % Wnf
        for o=1:n
            dWff_Wfn=abs(Wff(:,o)-Wfn)./Wfn;
            Rw_ff_fn(o)=wf_1(min(find(dWff_Wfn<=min(dWff_Wfn))));
            dWff_Wnf=abs(Wff(:,o)-Wnf(o))./Wnf(o);
            Rw_ff_nf(o)=wf_1(min(find(dWff_Wnf<=min(dWff_Wnf))));
            dWif_Wnf=abs(Wif(:,o)-Wnf(o))./Wnf(o);
            Rw_if_nf(o)=wi_1(min(find(dWif_Wnf<=min(dWif_Wnf))));
        end
        
        for i=1:n
            
            intmaxWin_Wnf(i)=sum((wi_1>Rw_in_nf(i)).*Fi_1.*derWin'.*wht_i_1);
            intmaxWni_Wnf(i)=sum((wi_2>Rw_ni_nf(i)).*Fi_2.*derWni'.*wht_i_2);
            intmaxWif_Wnf(i)=sum((wi_1>Rw_if_nf(i)).*Fi_1.*derWif(:,i).*wht_i_1);
            %intmaxWff_Wfn_Wnf(i)=sum((wf_1>Rw_ff_nf(i)).*Ff_1.*derWff(:,i).*wht_f_1);
            intmaxWff_Wfn_Wnf(i)=max(0,sum((wf_1<Rw_ff_fn(i)).*Ff_1.*derWff(:,i).*wht_f_1+(wf_1>=Rw_ff_fn(i)).*Ff_1.*derWfn'.*wht_f_1));
            
        end
        
        % Wni
        for o=1:n
            dWfi_Wfn=abs(Wfi(:,o)-Wfn)./Wfn;
            Rw_fi_fn(o)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        end
        
        for i=1:n
            
            intmaxWin_Wni(i)=sum((wi_1>Rw_in_ni(i)).*Fi_1.*derWin'.*wht_i_1);
            intmaxWnf_Wni(i)=sum((wf_2>Rw_nf_ni(i)).*Ff_2.*derWnf'.*wht_f_2);
            intmaxWii_Wni(i)=sum((wi_1>Rw_ii_ni(i)).*Fi_1.*derWii(:,i).*wht_i_1);
            %intmaxWfi_Wfn_Wni(i)=sum((wf_1>Rw_fi_ni(i)).*Ff_1.*derWfi(:,i).*wht_f_1);
            intmaxWfi_Wfn_Wni(i)=max(0,sum((wf_1<Rw_fi_fn(i)).*Ff_1.*derWfi(:,i).*wht_f_1+(wf_1>=Rw_fi_fn(i)).*Ff_1.*derWfn'.*wht_f_1));
            
        end
        
        
        for i=1:n
            for j=1:n
                
                % Wff
                
                intmaxWif_Wff(i,j)=sum((wi_1>Rw_if_ff(i,j)).*Fi_1.*derWif(:,j).*wht_i_1);
                intmaxWfi_Wff(i,j)=sum((wi_2>Rw_fi_ff(i,j)).*Fi_2.*derWfi(i,:)'.*wht_i_2);
                
            end
        end
        
        % Wii
        for o=1:n
            dWfi_Wfn=abs(Wfi(:,o)-Wfn)./Wfn;
            Rw_fi_fn(o)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
            dWif_Wnf=abs(Wif(o,:)'-Wnf)./Wnf;
            Rw_if_nf(o)=wf_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
        end
        
        for i=1:n
            for j=1:n
                
                % intmaxWfi_Wfn_Wii(i,j)=sum((wf_1>Rw_fi_ii(i,j)).*Ff_1.*derWfi(:,j).*wht_f_1);
                % intmaxWif_Wnf_Wii(i,j)=sum((wf_2>Rw_if_ii(i,j)).*Ff_2.*derWif(i,:)'.*wht_f_2);
                
                intmaxWfi_Wfn_Wii(i,j)=max(0,sum((wf_1<Rw_fi_fn(j)).*Ff_1.*derWfi(:,j).*wht_f_1+(wf_1>=Rw_fi_fn(j)).*Ff_1.*derWfn'.*wht_f_1));
                intmaxWif_Wnf_Wii(i,j)=max(0,sum((wf_2<Rw_if_nf(i)).*Ff_2.*derWif(i,:)'.*wht_f_2+(wf_2>=Rw_if_nf(i)).*Ff_2.*derWnf'.*wht_f_2));
           
              
            end
        end
        
        
        % Wfi
        for o=1:n
            dWff_Wnf=abs(Wff(o,:)'-Wnf)./Wnf;
            Rw_ff_nf(o)=wf_2(min(find(dWff_Wnf<=min(dWff_Wnf))));
        end
        
        for i=1:n
            for j=1:n
                intmaxWii_Wfi(i,j)=sum((wi_1>Rw_ii_fi(i,j)).*Fi_1.*derWii(:,j).*wht_i_1);
                % intmaxWff_Wnf_Wfi(i,j)=sum((wf_2>Rw_ff_fi(i,j)).*Ff_2.*derWff(i,:)'.*wht_f_2);
                intmaxWff_Wnf_Wfi(i,j)=max(0,sum((wf_2<Rw_ff_nf(i)).*Ff_2.*derWff(i,:)'.*wht_f_2+(wf_2>=Rw_ff_nf(i)).*Ff_2.*derWnf'.*wht_f_2));
                
            end
        end
        
        % Wif
        for o=1:n
            dWff_Wfn=abs(Wff(:,o)-Wfn)./Wfn;
            Rw_ff_fn(o)=wf_1(min(find(dWff_Wfn<=min(dWff_Wfn))));
        end
        
        for i=1:n
            for j=1:n
                
                intmaxWii_Wif(i,j)=sum((wi_2>Rw_ii_if(i,j)).*Fi_2.*derWii(i,:)'.*wht_i_2);
                % intmaxWff_Wfn_Wif(i,j)=sum((wf_1>Rw_ff_if(i,j)).*Ff_1.*derWff(:,j).*wht_f_1);
                intmaxWff_Wfn_Wif(i,j)=max(0,sum((wf_1<Rw_ff_fn(j)).*Ff_1.*derWff(:,j).*wht_f_1+(wf_1>=Rw_ff_fn(j)).*Ff_1.*derWfn'.*wht_f_1));
                
            end
        end
        
        % Wnn
        
        intmaxWfn_Wnn=sum((wf_1>Rw_fn_nn).*Ff_1.*derWfn'.*wht_f_1);
        intmaxWin_Wnn=sum((wi_1>Rw_in_nn).*Fi_1.*derWin'.*wht_i_1);
        intmaxWnf_Wnn=sum((wf_2>Rw_nf_nn).*Ff_2.*derWnf'.*wht_f_2);
        intmaxWni_Wnn=sum((wi_2>Rw_ni_nn).*Fi_2.*derWni'.*wht_i_2);
        
        
        B=di_2*(1-q_1)*(Wnn-Wni(1)) ...
            +di_2*q_1*intmaxWin_Wni(1) ...
            +lif_2*intmaxWnf_Wni(1) ...
            +lnf_1*intmaxWfi_Wfn_Wni(1)...
            +lni_1*intmaxWii_Wni(1);
        
        C=di_1*(1-q_2)*(Wnn-Win(1)) ...
            +di_1*q_2*intmaxWni_Win(1) ...
            +lif_1*intmaxWfn_Win(1) ...
            +lnf_2*intmaxWif_Wnf_Win(1) ...
            +lni_2*intmaxWii_Win(1);
                   
            % fixed across all simulations
%         b2=-(log(1-theta*(r*Wii(1,1)-gamma-C)))/theta-wi_1(1); % Wii(1,1)=Win(1)
%         b1=-(log(1-theta*(r*Wnn-gamma-B)))/theta-wi_2(1); % Wni(1)=Wnn    
%    


        %%%%% Value functions
        % solve in terms of F, reservation wages and lambdas
        
        % Wfn(w1) and Win(w1)
        for i=1:n
            
            eWfn(i)=(1/(r+df_1*(1-p_2)))*((1-exp(-theta*(wf_1(i)+b2)))/theta+a+df_1*(1-p_2)*Wnn ...
                +df_1*p_2*intmaxWni_Wfn(i) ...
                +lfi_1*intmaxWin_Wfn(i) ...
                +lnf_2*intmaxWff_Wnf_Wfn(i) ...
                +lni_2*intmaxWfi_Wfn(i))-Wfn(i);
            
            eWin(i)=(1/(r+di_1*(1-q_2)))*((1-exp(-theta*(wi_1(i)+b2)))/theta+gamma+di_1*(1-q_2)*Wnn ...
                +di_1*q_2*intmaxWni_Win(i) ...
                +lif_1*intmaxWfn_Win(i) ...
                +lnf_2*intmaxWif_Wnf_Win(i) ...
                +lni_2*intmaxWii_Win(i))-Win(i);

            
        end
        % Wnf(w2) and Wni(w2)
        for i=1:n
            
            eWnf(i)=(1/(r+df_2*(1-p_1)))*((1-exp(-theta*(wf_2(i)+b1)))/theta+a+df_2*(1-p_1)*Wnn ...
                +df_2*p_1*intmaxWin_Wnf(i) ...
                +lfi_2*intmaxWni_Wnf(i) ...
                +lnf_1*intmaxWff_Wfn_Wnf(i) ...
                +lni_1*intmaxWif_Wnf(i))-Wnf(i);
            
            eWni(i)=(1/(r+di_2*(1-q_1)))*((1-exp(-theta*(wi_2(i)+b1)))/theta+gamma+di_2*(1-q_1)*Wnn ...
                +di_2*q_1*intmaxWin_Wni(i) ...
                +lif_2*intmaxWnf_Wni(i) ...
                +lnf_1*intmaxWfi_Wfn_Wni(i) ...
                +lni_1*intmaxWii_Wni(i))-Wni(i);
            
        end
        
        for i=1:n
            for j=1:n
                % Wff(w1,w2)
                
                eWff(i,j)=(1/(r+df_1+df_2))*((1-exp(-theta*(wf_1(i)+wf_2(j))))/theta+a+df_1*Wnf(j)+df_2*Wfn(i) ...
                    +lfi_1*intmaxWif_Wff(i,j)...
                    +lfi_2*intmaxWfi_Wff(i,j))-Wff(i,j);
                
                % Wii(w1,w2)
                
                eWii(i,j)=(1/(r+di_1+di_2))*((1-exp(- theta*(wi_1(i)+wi_2(j))))/theta+gamma+di_1*Wni(j)+di_2*Win(i) ...
                    +lif_1*intmaxWfi_Wfn_Wii(i,j) ...
                    +lif_2*intmaxWif_Wnf_Wii(i,j))-Wii(i,j);
                
                
                % Wfi(w1,w2)
                
                eWfi(i,j)=(1/(r+df_1+di_2))*((1-exp(-theta*(wf_1(i)+wi_2(j))))/theta+a+df_1*Wni(j)+di_2*Wfn(i) ...
                    +lfi_1*intmaxWii_Wfi(i,j) ...
                    +lif_2*intmaxWff_Wnf_Wfi(i,j))-Wfi(i,j);
                
                
                % Wif(w1,w2)
                
                eWif(i,j)=(1/(r+di_1+df_2))*((1-exp(-theta*(wi_1(i)+wf_2(j))))/theta+a+di_1*Wnf(j)+df_2*Win(i) ...
                    +lif_1*intmaxWff_Wfn_Wif(i,j) ...
                    +lfi_2*intmaxWii_Wif(i,j))-Wif(i,j);
                
            end
        end
        % Wnn
        
        eWnn=(1/r)*((1-exp(-theta*(b1+b2)))/theta+gamma ...
            +lnf_1*intmaxWfn_Wnn ...
            +lni_1*intmaxWin_Wnn ...
            +lnf_2*intmaxWnf_Wnn ...
            +lni_2*intmaxWni_Wnn)-Wnn;

        
        % Update value functions
        
        Wfn=Wfn+omegaW*eWfn';
        Win=Win+omegaW*eWin';
        Wnf=Wnf+omegaW*eWnf';
        Wni=Wni+omegaW*eWni';
        Wff=Wff+omegaW*eWff;
        Wii=Wii+omegaW*eWii;
        Wfi=Wfi+omegaW*eWfi;
        Wif=Wif+omegaW*eWif;
        Wnn=Wni(1);eWnn=0;
        Wii(1,1)=Win(1);eWii(1,1)=0;
        
        critVF=[abs(eWfn'./Wfn);abs(eWin'./Win);abs(eWnf'./Wnf);abs(eWni'./Wni);max(abs(eWff./Wff))';max(abs(eWii./Wii))';max(abs(eWfi./Wfi))';max(abs(eWif./Wif))';abs(eWnn/Wnn)];
        %mean(critVF)
        iterVF=iterVF+1;
     
    end
    iterVF

end
